Sea ice thickness¶
In this notebook, we’ll show more detailed analysis and visualization of ICESat2 sea ice thickness measurements and the PIOMAS (Pan-Arctic Ice Ocean Modeling and Assimilation System) modelled sea ice thickness data. We’ll also highlight interactive mapping using the hvplot and holoviews package.
import xarray as xr # For working with gridded climate data
import pandas as pd
import numpy as np
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.misc_utils import restrictRegionally # Region restriction
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting
# Plotting dependencies
import cartopy.crs as ccrs
import holoviews as hv
from textwrap import wrap
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook
# Remove warnings to improve display
import warnings
warnings.filterwarnings('ignore')
1) Read in the data¶
We’ve included a function in the read_data_utils module that reads in the book dataset from the public google storage bucketas an xarray.Dataset object.
print(read_book_data.__doc__) # Print docstring
Read in data for ICESat2 jupyter book.
If the file does not already exist on the user's local drive, it is downloaded from the books google storage bucket (https://console.cloud.google.com/storage/browser/is2-pso-seaice)
The netcdf file is then read in as an xr.Dataset object
Args:
None
Returns:
book_ds (xr.Dataset): data
book_ds = read_book_data() # Read/download the data
print(book_ds) # Show dataset
<xarray.Dataset>
Dimensions: (time: 30, y: 448, x: 304)
Coordinates:
* time (time) datetime64[ns] 2018-11-01 ... 2021-04-01
longitude (y, x) float32 ...
latitude (y, x) float32 ...
region_mask (y, x) uint8 ...
Dimensions without coordinates: y, x
Data variables: (12/18)
freeboard (time, y, x) float32 ...
ice_density (time, y, x) float32 ...
ice_thickness (time, y, x) float32 ...
ice_thickness_smoothed (time, y, x) float64 ...
ice_thickness_unc (time, y, x) float32 ...
ice_type (time, y, x) float32 ...
... ...
piomas_ice_thickness (time, y, x) float32 ...
t2m (time, y, x) float32 ...
msdwlwrf (time, y, x) float32 ...
drifts_uT (time, y, x) float32 ...
drifts_vT (time, y, x) float32 ...
drifts_magnitude (time, y, x) float32 ...
Attributes:
description: Data used in ICESat-2 jupyter book
note: See individual data variables for references
creation date: 2021-09-20
Restrict data to a given region¶
We’ve built a region mask for the Arctic into the dataset used for this book; see the data wrangling notebook for more information. The region mask is included as a coordinate in the dataset, allowing us to easily access the different regions. You can restrict the data by selecting keys corresponding to regions of interest.
We are often interested in a region called the Inner Arctic, which is defined as the combined area of the Central Arctic, Beaufort Sea, Chukchi Sea, E Siberian Sea and the Laptev Sea. Figure from Petty et al., (2020):
Each region is defined by a specific integer key value. We can select the region/s defined by the key/s by using the restrictRegionally function, defined in the misc_utils module. Here, we’ll restrict the data to the Inner Arctic by selected the key values [10,11,12,13,15], corresponding to the five regions described above.
regions = pd.DataFrame({"labels":book_ds.region_mask.attrs["labels"]}, index=book_ds.region_mask.attrs["keys"])
print(regions)
labels
0 Lakes, extended coast
1 non-region oceans
2 Sea of Okhotsk and Japan
3 Bering Sea
4 Hudson Bay
5 Gulf of St. Lawrence
6 Baffin Bay, Davis Strait & Labrador Sea
7 Greenland Sea
8 Barents Seas
9 Kara Sea
10 Laptev Sea
11 East Siberian Sea
12 Chukchi Sea
13 Beaufort Sea
14 Canadian Archipelago
15 Arctic Ocean
20 Land
21 Coast
The regions defining the Inner Arctic are associated with the key values 10, 11, 12, 13, and 15. Here, we’ll use the restrictRegionally function defined in the misc_utils module to mask out data outside of that region.
innerArcticKeys = [10,11,12,13,15] #Inner Arctic
book_ds_innerArctic = restrictRegionally(book_ds, innerArcticKeys)
Regions selected: Inner Arctic
On the map below, you can see how the data we’ve called the function on has been cut to show only gridcells in the Inner Arctic!
date = "Feb 2020"
fig, axes = plt.subplots(1, 2, figsize=(5,5), subplot_kw={'projection':ccrs.NorthPolarStereo(central_longitude=-45)})
im1 = book_ds["ice_thickness_smoothed"].sel(time=date)[0].plot.pcolormesh(ax=axes[0], x='longitude', y='latitude', vmin=0, vmax=4, extend='both', transform=ccrs.PlateCarree(), add_colorbar=False)
im2 = book_ds_innerArctic["ice_thickness_smoothed"].sel(time=date)[0].plot.pcolormesh(ax=axes[1], x='longitude', y='latitude', vmin=0, vmax=4, extend='both', transform=ccrs.PlateCarree(), add_colorbar=False)
plt.colorbar(im1, cax=fig.add_axes([0.93, 0.25, 0.025, 0.5]), extend='both', label="sea ice thickness (m)")
axes[0].set_title("\n".join(wrap("Sea ice thickness before restricting, "+date, 32)), fontsize=8)
axes[1].set_title("\n".join(wrap("Inner Arctic sea ice thickness, "+date, 32)), fontsize=8)
axes[0].coastlines()
axes[1].coastlines()
plt.show()
Set years over which to perform the analysis¶
Setting years = [2018,2019,2020] indicates that you want to perform the analyses in the notebook for three winter season: winter 2018-19, 2019-20, and 2020-21
years = [2018,2019,2020]
3) ICESat-2 interpolated sea ice thickness maps¶
Here, we’ll overlay sea ice thickness data over a map of the Arctic. To generate the map, we’ll use interactive plotting functions based off the hvplot python package. To see more information on the functions used, click on the plus buttons to see the function documentation, or go into the plotting_utils module to see the code.
print(interactiveArcticMaps.__doc__) # Print docstring
Generative one or more interactive maps
Using the argument "slide", the user can set whether each map should be displayed together, or displayed in the form of a slider
To show each map together (no slider), set slider=False
Args:
da (xr.Dataset or xr.DataArray): data to restrict by time; must contain "time" as a coordinate
clabel (str, optional): colorbar label (default to "long_name" and "units" if given in attributes of da)
cmap (str, optional): matplotlib colormap to use (default to "viridis")
colorbar (bool, optional): show colorbar? (default to True)
vmin (float, optional): minimum on colorbar (default to 1st percentile)
vmax (float, optional): maximum on colorbar (default to 99th percentile)
title (str, optional): main title to give plot (default to no title)
ylim (tuple, optional): limits of yaxis in the form min latitude, max latitude (default to (60,90))
frame_width (int, optional): width of frame. sets figure size of each map (default to 250)
slider (bool, optional): if da has more than one time coordinate, display maps with a slider? (default to True)
cols (int, optional): how many columns to show before wrapping, if da has more than one time coordinate (default to 3)
Returns:
pl (Holoviews map)
pl = interactiveArcticMaps(book_ds["ice_thickness_smoothed"], colorbar=True, slider=True, frame_width=500)
pl
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
/var/folders/8v/cr06mz0n3bjd5jr12_9v6n200000gn/T/ipykernel_35719/2060953563.py in <module>
1 pl = interactiveArcticMaps(book_ds["ice_thickness_smoothed"], colorbar=True, slider=True, frame_width=500)
----> 2 pl
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/displayhook.py in __call__(self, result)
260 self.start_displayhook()
261 self.write_output_prompt()
--> 262 format_dict, md_dict = self.compute_format_data(result)
263 self.update_user_ns(result)
264 self.fill_exec_result(result)
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/displayhook.py in compute_format_data(self, result)
149
150 """
--> 151 return self.shell.display_formatter.format(result)
152
153 # This can be set to True by the write_output_prompt method in a subclass
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
148 return {}, {}
149
--> 150 format_dict, md_dict = self.mimebundle_formatter(obj, include=include, exclude=exclude)
151
152 if format_dict or md_dict:
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/decorator.py in fun(*args, **kw)
230 if not kwsyntax:
231 args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)
233 fun.__name__ = func.__name__
234 fun.__doc__ = func.__doc__
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
222 """show traceback on failed format call"""
223 try:
--> 224 r = method(self, *args, **kwargs)
225 except NotImplementedError:
226 # don't warn on NotImplementedErrors
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/IPython/core/formatters.py in __call__(self, obj, include, exclude)
968
969 if method is not None:
--> 970 return method(include=include, exclude=exclude)
971 return None
972 else:
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/core/dimension.py in _repr_mimebundle_(self, include, exclude)
1315 combined and returned.
1316 """
-> 1317 return Store.render(self)
1318
1319
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/core/options.py in render(cls, obj)
1403 data, metadata = {}, {}
1404 for hook in hooks:
-> 1405 ret = hook(obj)
1406 if ret is None:
1407 continue
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in pprint_display(obj)
280 if not ip.display_formatter.formatters['text/plain'].pprint:
281 return None
--> 282 return display(obj, raw_output=True)
283
284
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in display(obj, raw_output, **kwargs)
256 elif isinstance(obj, (HoloMap, DynamicMap)):
257 with option_state(obj):
--> 258 output = map_display(obj)
259 elif isinstance(obj, Plot):
260 output = render(obj)
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in wrapped(element)
144 try:
145 max_frames = OutputSettings.options['max_frames']
--> 146 mimebundle = fn(element, max_frames=max_frames)
147 if mimebundle is None:
148 return {}, {}
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in map_display(vmap, max_frames)
204 return None
205
--> 206 return render(vmap)
207
208
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/ipython/display_hooks.py in render(obj, **kwargs)
66 renderer = renderer.instance(fig='png')
67
---> 68 return renderer.components(obj, **kwargs)
69
70
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/renderer.py in components(self, obj, fmt, comm, **kwargs)
408 doc = Document()
409 with config.set(embed=embed):
--> 410 model = plot.layout._render_model(doc, comm)
411 if embed:
412 return render_model(model, comm)
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/viewable.py in _render_model(self, doc, comm)
453 if comm is None:
454 comm = state._comm_manager.get_server_comm()
--> 455 model = self.get_root(doc, comm)
456
457 if config.embed:
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/viewable.py in get_root(self, doc, comm, preprocess)
510 """
511 doc = init_doc(doc)
--> 512 root = self._get_model(doc, comm=comm)
513 if preprocess:
514 self._preprocess(root)
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/layout/base.py in _get_model(self, doc, root, parent, comm)
120 if root is None:
121 root = model
--> 122 objects = self._get_objects(model, [], doc, root, comm)
123 props = dict(self._init_params(), objects=objects)
124 model.update(**self._process_param_change(props))
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/layout/base.py in _get_objects(self, model, old_objects, doc, root, comm)
110 else:
111 try:
--> 112 child = pane._get_model(doc, root, model, comm)
113 except RerenderError:
114 return self._get_objects(model, current_objects[:i], doc, root, comm)
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/pane/holoviews.py in _get_model(self, doc, root, parent, comm)
261 if p in Layoutable.param and p != 'name'}
262 child_pane = self._panes.get(backend, Pane)(state, **kwargs)
--> 263 self._update_plot(plot, child_pane)
264 model = child_pane._get_model(doc, root, parent, comm)
265 if ref in self._plots:
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/panel/pane/holoviews.py in _update_plot(self, plot, pane)
197 if plot.comm or state._unblocked(plot.document):
198 with unlocked():
--> 199 plot.update(key)
200 if plot.comm and 'embedded' not in plot.root.tags:
201 plot.push()
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/plot.py in update(self, key)
981 if len(self) == 1 and ((key == 0) or (key == self.keys[0])) and not self.drawn:
982 return self.initialize_plot()
--> 983 item = self.__getitem__(key)
984 self.traverse(lambda x: setattr(x, '_updated', True))
985 return item
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/plot.py in __getitem__(self, frame)
444 if not isinstance(frame, tuple):
445 frame = self.keys[frame]
--> 446 self.update_frame(frame)
447 return self.state
448
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py in update_frame(self, key, ranges, element)
2452 not subplot in self.dynamic_subplots):
2453 continue
-> 2454 subplot.update_frame(key, ranges, element=el)
2455
2456 if not self.batched and isinstance(self.hmap, DynamicMap) and items:
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py in update_frame(self, key, ranges, plot, element)
1541 self._update_hover(element)
1542
-> 1543 self._update_glyphs(element, ranges, self.style[self.cyclic_index])
1544 self._execute_hooks(element)
1545
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/element.py in _update_glyphs(self, element, ranges, style)
1448 data, mapping, style = self.get_batched_data(element, ranges)
1449 else:
-> 1450 data, mapping, style = self.get_data(element, ranges, style)
1451
1452 # Include old data if source static
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/geoviews/plotting/bokeh/plot.py in get_data(self, element, ranges, style)
171 if self._project_operation and self.geographic:
172 element = self._project_operation(element, projection=self.projection)
--> 173 return super(GeoPlot, self).get_data(element, ranges, style)
174
175
~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/holoviews/plotting/bokeh/raster.py in get_data(self, element, ranges, style)
270 yc.append(ys.mean())
271 else:
--> 272 mask.append(False)
273
274 data = {'xs': XS, 'ys': YS, z.name: zvals[np.array(mask)]}
KeyboardInterrupt:
print(interactive_winter_mean_maps.__doc__) # Print docstring
Generate interactive maps of winter mean data
Note: this function builds off the functions get_winter_data and interactiveArcticMaps.
Args:
da (xr.Dataset or xr.DataArray): data to restrict by time; must contain "time" as a coordinate
years (list of str): years over which to compute mean (default to unique years in the dataset)
start_month (str, optional): first month in winter (default to September)
end_month (str, optional): second month in winter; this is the following calender year after start_month (default to April)
force_complete_season (bool, optional): require that winter season returns data if and only if all months have data? i.e. if Sep and Oct have no data, return nothing even if Nov-Apr have data? (default to False)
clabel (str, optional): colorbar label (default to "long_name" and "units" if given in attributes of da)
cmap (str, optional): matplotlib colormap to use (default to "viridis")
colorbar (bool, optional): show colorbar? (default to True)
vmin (float, optional): minimum on colorbar (default to 0)
vmax (float, optional): maximum on colorbar (default to 4)
title (str, optional): main title to give plot (default to no title)
ylim (tuple, optional): limits of yaxis in the form min latitude, max latitude (default to (60,90))
frame_width (int, optional): width of frame. sets figure size of each map (default to 250)
slider (bool, optional): if da has more than one time coordinate, display maps with a slider? (default to True)
cols (int, optional): how many columns to show before wrapping, if da has more than one time coordinate (default to 3)
Returns:
pl_means (Holoviews map)
interactive_winter_mean_maps(book_ds_innerArctic["ice_thickness_smoothed"],
years=years,
slider=True,
frame_width=500,
title="ICESat-2 sea ice thickness winter means")
4) Differences in winter sea ice thickness¶
Here, we compare the sea ice thickness across winter seasons and generate a map of the gridcells differences to identidy any regions of particularly high differences in either direction. This allows us to get an idea of how and where two winter seasons differed in thickness.
Because ICESat-2 didn’t start collecting data until November 2018, we’ll show means for November through April, instead of Sep through April, which we would normally define as the winter season.
The gridcell difference map shows the middle plot minus the first plot, so a positive (colored red in the map) difference would indicate that winter 2019-2020 had a higher sea ice thickness value than winter 2018-2019, and a negative (blue) difference would indicate a lower sea ice thickness value.
is2_winter_means = compute_gridcell_winter_means(book_ds_innerArctic["ice_thickness_smoothed"], years=years, start_month="Nov")
pio_winter_means = compute_gridcell_winter_means(book_ds_innerArctic["piomas_ice_thickness"], years=years, start_month="Nov")
for i in range(len(years)):
data1 = is2_winter_means.isel(time=i)
data2 = pio_winter_means.isel(time=i)
pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=data1.time.values.item())
pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title=data2.time.values.item())
diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference")
data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
display(data_pl)
print("\n")